Bayesian GLM Part7

Author

Murray Logan

Published

08/07/2025

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(DHARMa) # for residual diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(tidybayes) # for more tidying outputs
library(HDInterval) # for HPD intervals
library(ggeffects) # for partial plots
library(broom.mixed) # for summarising models
library(posterior) # for posterior draws
library(ggeffects) # for partial effects plots
library(patchwork) # for multi-panel figures
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(geoR) # framework for stats, modelling and visualisation
library(modelsummary) # for data and model summaries
theme_set(theme_grey()) # put the default ggplot theme back
source("helperFunctions.R")

2 Scenario

Mitchell et al. (2019) presented a data set in which male guppy fish where placed individually into long, narrow and shallow tanks that predominately only permitted the fish to swim in two dimensions. The activity rate (cumulative distance swam) of each guppy fish was remotely recorded via tracking software during 20 minute trials conducted each day for 14 days. The order in which individual guppy fish were recorded each day was randomised and the time of day of each trial was recorded.

The researchers were interested in whether the the guppy fish become habituated to the experimental conditions and altered their activity patterns over time.

Figure 1: guppy fish

The data are in the file mitchell.csv in the data folder.

Table 1: Format of the mitchell.csv data file
ID TOD Day Dist_moved
A39 0.39 7 3439.93
A3 0.39 7 2795.08
A13 0.39 7 5258.79
... ...
Table 2: Description of the variables in the mitchell data file
ID: Individual guppy fish ID - Random variable
TOD: Time of day of observation as a proportion - Predictor variable
Day: Day post entry into tank - Predictor variable
Dist_moved: Total distance (cm) moved in 20 min trial - Response variable

3 Read in the data

mitchell <- read_csv("../data/mitchell.csv", trim_ws = TRUE)
Rows: 1626 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): ID, Date, StartDate
dbl  (5): Mass, Day, TOD, Dist_moved, Dist
time (1): TimeOfDay

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(mitchell)
Rows: 1,626
Columns: 9
$ ID         <chr> "A39", "A3", "A13", "A23", "A59", "A69", "A41", "A9", "A35"…
$ TimeOfDay  <time> 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32, 09:28:32…
$ Date       <chr> "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2018", "6/03/2…
$ StartDate  <chr> "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb-18", "27-Feb…
$ Mass       <dbl> 0.066, 0.072, 0.105, 0.080, 0.077, 0.092, 0.073, 0.075, 0.0…
$ Day        <dbl> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
$ TOD        <dbl> 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.39, 0.42, 0.42, 0.42,…
$ Dist_moved <dbl> 3439.930, 2795.080, 5258.790, 3402.740, 7672.660, 2400.660,…
$ Dist       <dbl> 58.65092, 52.86852, 72.51752, 58.33301, 87.59372, 48.99653,…
## Explore the first 6 rows of the data
head(mitchell)
# A tibble: 6 × 9
  ID    TimeOfDay Date      StartDate  Mass   Day   TOD Dist_moved  Dist
  <chr> <time>    <chr>     <chr>     <dbl> <dbl> <dbl>      <dbl> <dbl>
1 A39   09:28:32  6/03/2018 27-Feb-18 0.066     7  0.39      3440.  58.7
2 A3    09:28:32  6/03/2018 27-Feb-18 0.072     7  0.39      2795.  52.9
3 A13   09:28:32  6/03/2018 27-Feb-18 0.105     7  0.39      5259.  72.5
4 A23   09:28:32  6/03/2018 27-Feb-18 0.08      7  0.39      3403.  58.3
5 A59   09:28:32  6/03/2018 27-Feb-18 0.077     7  0.39      7673.  87.6
6 A69   09:28:32  6/03/2018 27-Feb-18 0.092     7  0.39      2401.  49.0
str(mitchell)
spc_tbl_ [1,626 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ ID        : chr [1:1626] "A39" "A3" "A13" "A23" ...
 $ TimeOfDay : 'hms' num [1:1626] 09:28:32 09:28:32 09:28:32 09:28:32 ...
  ..- attr(*, "units")= chr "secs"
 $ Date      : chr [1:1626] "6/03/2018" "6/03/2018" "6/03/2018" "6/03/2018" ...
 $ StartDate : chr [1:1626] "27-Feb-18" "27-Feb-18" "27-Feb-18" "27-Feb-18" ...
 $ Mass      : num [1:1626] 0.066 0.072 0.105 0.08 0.077 0.092 0.073 0.075 0.09 0.077 ...
 $ Day       : num [1:1626] 7 7 7 7 7 7 7 7 7 7 ...
 $ TOD       : num [1:1626] 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.42 0.42 0.42 ...
 $ Dist_moved: num [1:1626] 3440 2795 5259 3403 7673 ...
 $ Dist      : num [1:1626] 58.7 52.9 72.5 58.3 87.6 ...
 - attr(*, "spec")=
  .. cols(
  ..   ID = col_character(),
  ..   TimeOfDay = col_time(format = ""),
  ..   Date = col_character(),
  ..   StartDate = col_character(),
  ..   Mass = col_double(),
  ..   Day = col_double(),
  ..   TOD = col_double(),
  ..   Dist_moved = col_double(),
  ..   Dist = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
mitchell |> datawizard::data_codebook()
mitchell (1626 rows and 9 variables, 9 shown)

ID | Name       | Type      | Missings |            Values |           N
---+------------+-----------+----------+-------------------+------------
1  | ID         | character | 0 (0.0%) |                A1 |  14 ( 0.9%)
   |            |           |          |               A11 |  14 ( 0.9%)
   |            |           |          |               A13 |  14 ( 0.9%)
   |            |           |          |               A15 |  14 ( 0.9%)
   |            |           |          |               A17 |  14 ( 0.9%)
   |            |           |          |               A19 |  14 ( 0.9%)
   |            |           |          |               A21 |  14 ( 0.9%)
   |            |           |          |               A23 |  14 ( 0.9%)
   |            |           |          |               A25 |  14 ( 0.9%)
   |            |           |          |               A27 |  14 ( 0.9%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
2  | TimeOfDay  | numeric   | 0 (0.0%) |          07:57:31 |   7 ( 0.4%)
   |            |           |          |          08:08:00 |   7 ( 0.4%)
   |            |           |          |          08:16:18 |   2 ( 0.1%)
   |            |           |          |          08:17:16 |   8 ( 0.5%)
   |            |           |          |          08:17:56 |   5 ( 0.3%)
   |            |           |          |          08:18:36 |   7 ( 0.4%)
   |            |           |          |          08:20:22 |   6 ( 0.4%)
   |            |           |          |          08:23:28 |   8 ( 0.5%)
   |            |           |          |          08:26:07 |   5 ( 0.3%)
   |            |           |          |          08:30:43 |  10 ( 0.6%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
3  | Date       | character | 0 (0.0%) |         1/04/2018 |  40 ( 2.5%)
   |            |           |          |         1/05/2018 |  40 ( 2.5%)
   |            |           |          |        10/03/2018 |  39 ( 2.4%)
   |            |           |          |        10/04/2018 |  40 ( 2.5%)
   |            |           |          |        11/03/2018 |  39 ( 2.4%)
   |            |           |          |        11/04/2018 |  40 ( 2.5%)
   |            |           |          |        12/03/2018 |  39 ( 2.4%)
   |            |           |          |        12/04/2018 |  40 ( 2.5%)
   |            |           |          |        13/03/2018 |  39 ( 2.4%)
   |            |           |          |        14/03/2018 |  39 ( 2.4%)
   |            |           |          |             (...) |            
---+------------+-----------+----------+-------------------+------------
4  | StartDate  | character | 0 (0.0%) |         17-Apr-18 | 560 (34.4%)
   |            |           |          |         23-Mar-18 | 520 (32.0%)
   |            |           |          |         27-Feb-18 | 546 (33.6%)
---+------------+-----------+----------+-------------------+------------
5  | Mass       | numeric   | 0 (0.0%) |      [0.05, 0.12] |        1626
---+------------+-----------+----------+-------------------+------------
6  | Day        | numeric   | 0 (0.0%) |           [7, 20] |        1626
---+------------+-----------+----------+-------------------+------------
7  | TOD        | numeric   | 0 (0.0%) |      [0.37, 0.56] |        1626
---+------------+-----------+----------+-------------------+------------
8  | Dist_moved | numeric   | 0 (0.0%) | [101.13, 11833.9] |        1626
---+------------+-----------+----------+-------------------+------------
9  | Dist       | numeric   | 0 (0.0%) |   [10.06, 108.78] |        1626
------------------------------------------------------------------------
mitchell |> modelsummary::datasummary_skim()
Warning: These variables were omitted because they include more than 50 levels:
ID.
Unique Missing Pct. Mean SD Min Median Max Histogram
Mass 51 0 0.1 0.0 0.0 0.1 0.1
Day 14 0 13.5 4.1 7.0 13.0 20.0
TOD 97 0 0.5 0.0 0.4 0.4 0.6
Dist_moved 1625 0 3917.0 2141.8 101.1 3811.0 11833.9
Dist 1625 0 59.9 18.1 10.1 61.7 108.8
N %
Date 1/04/2018 40 2.5
1/05/2018 40 2.5
10/03/2018 39 2.4
10/04/2018 40 2.5
11/03/2018 39 2.4
11/04/2018 40 2.5
12/03/2018 39 2.4
12/04/2018 40 2.5
13/03/2018 39 2.4
14/03/2018 39 2.4
15/03/2018 39 2.4
16/03/2018 39 2.4
17/03/2018 39 2.4
18/03/2018 39 2.4
19/03/2018 39 2.4
2/04/2018 40 2.5
2/05/2018 40 2.5
24/04/2018 40 2.5
25/04/2018 40 2.5
26/04/2018 40 2.5
27/04/2018 40 2.5
28/04/2018 40 2.5
29/04/2018 40 2.5
3/04/2018 40 2.5
3/05/2018 40 2.5
30/03/2018 40 2.5
30/04/2018 40 2.5
31/03/2018 40 2.5
4/04/2018 40 2.5
4/05/2018 40 2.5
5/04/2018 40 2.5
5/05/2018 40 2.5
6/03/2018 39 2.4
6/04/2018 40 2.5
6/05/2018 40 2.5
7/03/2018 39 2.4
7/05/2018 40 2.5
8/03/2018 39 2.4
8/04/2018 40 2.5
9/03/2018 39 2.4
9/04/2018 40 2.5
StartDate 17-Apr-18 560 34.4
23-Mar-18 520 32.0
27-Feb-18 546 33.6
mitchell |> modelsummary::datasummary_skim(by = "Day")
Warning: These variables were omitted because they include more than 50 levels:
ID.
Day Unique Missing Pct. Mean SD Min Median Max Histogram
Mass 7.0 51 0 0.1 0.0 0.0 0.1 0.1
8.0 51 0 0.1 0.0 0.0 0.1 0.1
9.0 51 0 0.1 0.0 0.0 0.1 0.1
10.0 51 0 0.1 0.0 0.0 0.1 0.1
11.0 51 0 0.1 0.0 0.0 0.1 0.1
12.0 51 0 0.1 0.0 0.0 0.1 0.1
13.0 51 0 0.1 0.0 0.0 0.1 0.1
14.0 51 0 0.1 0.0 0.0 0.1 0.1
15.0 45 0 0.1 0.0 0.0 0.1 0.1
16.0 51 0 0.1 0.0 0.0 0.1 0.1
17.0 51 0 0.1 0.0 0.0 0.1 0.1
18.0 51 0 0.1 0.0 0.0 0.1 0.1
19.0 51 0 0.1 0.0 0.0 0.1 0.1
20.0 51 0 0.1 0.0 0.0 0.1 0.1
TOD 7.0 12 0 0.5 0.0 0.4 0.5 0.5
8.0 11 0 0.5 0.1 0.4 0.4 0.6
9.0 18 0 0.5 0.1 0.4 0.4 0.6
10.0 18 0 0.5 0.0 0.4 0.5 0.5
11.0 19 0 0.4 0.0 0.4 0.4 0.5
12.0 19 0 0.4 0.0 0.4 0.4 0.5
13.0 20 0 0.5 0.1 0.4 0.4 0.5
14.0 21 0 0.5 0.1 0.4 0.4 0.5
15.0 13 0 0.5 0.1 0.4 0.5 0.6
16.0 17 0 0.5 0.0 0.4 0.4 0.5
17.0 17 0 0.5 0.0 0.4 0.5 0.6
18.0 20 0 0.5 0.1 0.4 0.5 0.5
19.0 17 0 0.5 0.0 0.4 0.4 0.6
20.0 18 0 0.5 0.0 0.4 0.5 0.5
Dist_moved 7.0 119 0 3766.9 2207.2 101.1 3402.7 9037.0
8.0 119 0 3609.6 2066.8 269.8 3589.3 8343.5
9.0 119 0 3468.2 1996.3 251.3 3396.0 8581.9
10.0 119 0 3525.7 1908.6 472.4 3256.8 8343.1
11.0 119 0 3603.1 2032.4 341.6 3389.6 8508.7
12.0 119 0 3682.0 2081.9 577.9 3425.5 9167.5
13.0 119 0 3918.3 2181.0 415.2 3691.7 8840.6
14.0 119 0 4018.1 2107.9 566.2 3928.4 8881.7
15.0 79 0 4031.3 2212.8 685.7 3851.9 9619.0
16.0 119 0 4199.5 2241.9 438.1 4215.6 10238.8
17.0 119 0 4234.3 2222.7 549.3 4107.0 9439.1
18.0 119 0 4315.3 2128.9 563.4 4684.7 9420.1
19.0 119 0 4238.4 2202.5 656.7 4053.9 9186.6
20.0 119 0 4266.2 2220.8 433.4 4454.1 11833.9
Dist 7.0 119 0 58.3 19.4 10.1 58.3 95.1
8.0 119 0 57.3 18.1 16.4 59.9 91.3
9.0 119 0 56.2 17.7 15.9 58.3 92.6
10.0 119 0 56.9 16.9 21.7 57.1 91.3
11.0 119 0 57.3 18.0 18.5 58.2 92.2
12.0 119 0 58.1 17.6 24.0 58.5 95.7
13.0 119 0 59.8 18.6 20.4 60.8 94.0
14.0 119 0 60.9 17.7 23.8 62.7 94.2
15.0 79 0 60.8 18.4 26.2 62.1 98.1
16.0 119 0 62.2 18.4 20.9 64.9 101.2
17.0 119 0 62.6 18.0 23.4 64.1 97.2
18.0 119 0 63.3 17.6 23.7 68.4 97.1
19.0 119 0 62.6 18.0 25.6 63.7 95.8
20.0 119 0 62.8 18.1 20.8 66.7 108.8
N %
Date 1/04/2018 40 2.5
1/05/2018 40 2.5
10/03/2018 39 2.4
10/04/2018 40 2.5
11/03/2018 39 2.4
11/04/2018 40 2.5
12/03/2018 39 2.4
12/04/2018 40 2.5
13/03/2018 39 2.4
14/03/2018 39 2.4
15/03/2018 39 2.4
16/03/2018 39 2.4
17/03/2018 39 2.4
18/03/2018 39 2.4
19/03/2018 39 2.4
2/04/2018 40 2.5
2/05/2018 40 2.5
24/04/2018 40 2.5
25/04/2018 40 2.5
26/04/2018 40 2.5
27/04/2018 40 2.5
28/04/2018 40 2.5
29/04/2018 40 2.5
3/04/2018 40 2.5
3/05/2018 40 2.5
30/03/2018 40 2.5
30/04/2018 40 2.5
31/03/2018 40 2.5
4/04/2018 40 2.5
4/05/2018 40 2.5
5/04/2018 40 2.5
5/05/2018 40 2.5
6/03/2018 39 2.4
6/04/2018 40 2.5
6/05/2018 40 2.5
7/03/2018 39 2.4
7/05/2018 40 2.5
8/03/2018 39 2.4
8/04/2018 40 2.5
9/03/2018 39 2.4
9/04/2018 40 2.5
StartDate 17-Apr-18 560 34.4
23-Mar-18 520 32.0
27-Feb-18 546 33.6

4 Data preparation

Let start by declaring the categorical variables and random effect as factors.

mitchell <- mitchell |>
  mutate(
    ID = factor(ID),
    fDay = factor(Day)
  )

5 Exploratory data analysis

As the response represents a strictly positive real number (total distance moved), it would make sense to start by considering a Gamma model.

Note, also that the response variable are relatively large distances (as the units of distance used were centimeters). It might be worth dividing these values by 100 or even 1000 during the modelling process.

Model formula: \[ y_i/1000 \sim{} \mathcal{Gamma}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of day and time of day on the total distance moved during the trial. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual guppies.

mitchell |> ggplot(aes(y = Dist_moved)) +
  geom_boxplot()

## mitchell |> ggplot(aes(y = Dist)) +
##   geom_boxplot()
mitchell |> ggplot(aes(y = Dist_moved)) +
  geom_boxplot() +
  scale_y_log10()

mitchell |> ggplot(aes(y = TOD)) +
  geom_boxplot() +
  scale_y_log10()

mitchell |> ggplot(aes(y = Day)) +
  geom_boxplot()

## mitchell |> ggplot(aes(y = Dist, x = Day)) +
##   geom_point() +
##   geom_line(aes(group = ID))
mitchell |> ggplot(aes(y = Dist_moved, x = TOD)) +
  geom_point() +
  geom_line(aes(group = ID))

mitchell |> ggplot(aes(y = Dist, x = Day)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  geom_smooth(method = "lm") +
  facet_wrap(~ID)
`geom_smooth()` using formula = 'y ~ x'

mitchell |> ggplot(aes(y = Dist, x = TOD)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  geom_smooth(method = "lm") +
  facet_wrap(~ID)
`geom_smooth()` using formula = 'y ~ x'

mitchell |> ggplot(aes(y = Day, x = TOD)) +
  geom_point()

6 Fit the model

mitchell_sub <- mitchell |>
  ## filter(ID == "A63") |>
  filter(ID == "C69") |>
  droplevels()

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

mitchell.form <- bf(I(Dist_moved/1000) ~ scale(Day) + scale(TOD),
                family=Gamma(link='log'))
options(width=150)
mitchell.form |> get_prior(data = mitchell_sub)
                  prior     class     coef group resp dpar nlpar lb ub       source
                 (flat)         b                                           default
                 (flat)         b scaleDay                             (vectorized)
                 (flat)         b scaleTOD                             (vectorized)
 student_t(3, 1.3, 2.5) Intercept                                           default
      gamma(0.01, 0.01)     shape                                 0         default
options(width=80)

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

  • \(\beta_0\): normal centred at 134 with a standard deviation of 65
    • mean of 134: since median(log(mitchell$Dist_moved/1000))
    • sd pf 0.7: since mad(log(mitchell$Dist_moved/1000))
  • \(\beta_1\): normal centred at 0 with a standard deviation of 2
    • sd of 1: since mad(log(Dist_moved/1000))/mad(scale(Day))
  • \(\sigma\): half-t centred at 0 with a standard deviation of 65 OR
    • sd pf 65: since mad(log(mitchell$Dist_moved/1000))
  • \(\sigma\): gamma with shape parameters of 2 and 1
  • \(\beta_0\): normal centred at 1.5 with a standard deviation of 1.5
    • mean of 1.5: since mean(log(mitchell$Dist_moved/1000)) or mean(asinh(mitchell$Dist_moved/(2*1))/log(exp(1)))
    • sd of 1.5: since sd(log(mitchell$Dist_moved/1000)) or sd(asinh(mitchell$Dist_moved/(2*1))/log(exp(1)))
mitchell_sub |>
  mutate(Dist = Dist) |>
  summarise(
    Int_mu = median(Dist),
    Int_sd = mad(Dist),
    b1_sd = mad(Dist) / mad(scale(Day)),
    b2_sd = mad(Dist) / mad(scale(TOD))
  )
# A tibble: 1 × 4
  Int_mu Int_sd b1_sd b2_sd
   <dbl>  <dbl> <dbl> <dbl>
1   61.5   13.4  10.8  9.97
## mitchell |>
##   mutate(lDist = log(Dist_moved/1)) |>
##   summarise(
##     Int_mu = median(lDist),
##     Int_sd = mad(lDist),
##     b1_sd = mad(lDist)/mad(scale(Day)),
##     b2_sd = mad(lDist)/mad(scale(TOD))
##     )
## standist::visualize("normal(1.4, 1)", xlim=c(0,20))
## standist::visualize("student_t(3, 0, 1)",
##                     xlim=c(-10,25))
## priors <- prior(normal(1.4, 1), class = 'Intercept') +
##     prior(normal(0, 1), class = 'b') +
##     prior(student_t(3, 0, 1), class = 'sd')
##     prior(gamma(2, 1), class = 'shape')

priors <- prior(normal(60, 25), class = "Intercept") +
  prior(normal(0, 15), class = "b") +
  ## prior(student_t(3, 0, 25), class = 'sd') +
  prior(student_t(3, 0, 25), class = "sigma")
## prior(gamma(2, 1), class = 'shape')
## mitchell.form <- bf(Dist_moved ~ scale(Day) + scale(TOD) + (scale(Day) + scale(TOD)|ID),
##                 ## family=Gamma(link='log'))
##                 family=gaussian())
mitchell.form <- bf(Dist ~ scale(Day) + scale(TOD),
  ## family=Gamma(link='log'))
  family = gaussian()
)
mitchell.brm2 <- brm(mitchell.form,
  data = mitchell_sub,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 100,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)
Compiling Stan program...
Start sampling
mitchell.brm2 |>
  conditional_effects("Day") |>
  plot(points = TRUE)

mitchell.brm2 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm2 |>
  conditional_effects("TOD") |>
  plot(points = TRUE)

mitchell.brm2 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.

mitchell.brm3 <- update(mitchell.brm2,
  sample_prior = "yes",
  chains = 3, cores = 3,
  refresh = 100
)
The desired updates require recompiling the model
Compiling Stan program...
Start sampling
mitchell.brm3 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm3 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm3 |> SUYR_prior_and_posterior()

7 MCMC sampling diagnostics

The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
mitchell.brm3$fit |> stan_trace()

mitchell.brm3$fit |> stan_trace(inc_warmup = TRUE)

The chains appear well mixed and very similar

  • stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
mitchell.brm3$fit |> stan_ac()

There is no evidence of auto-correlation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
mitchell.brm3$fit |> stan_rhat()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

mitchell.brm3$fit |> stan_ess()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

mitchell.brm3$fit |> stan_dens(separate_chains = TRUE)

8 Model validation

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()
bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped
  • dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
mitchell.brm3 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to represent the shape of the observed data reasonably well

  • error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.
mitchell.brm3 |> pp_check(type = "error_scatter_avg")
Using all posterior draws for ppc type 'error_scatter_avg' by default.

The predictive error seems to be related to the predictor - the model performs poorest at higher recruitments.

  • intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
mitchell.brm3 |> pp_check(type = "intervals")
Using all posterior draws for ppc type 'intervals' by default.

## mitchell.brm3 |> pp_check(group='DENSITY', type='intervals')

Whilst the modelled predictions do a reasonable job of representing the observed data, the observed data do appear to be more varied than the model is representing.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(mitchell.brm3)

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

  • simulated (predicted) responses associated with each observation.
  • observed values
  • fitted (predicted) responses (averaged) associated with each observation
preds <- mitchell.brm3 |> posterior_predict(nsamples = 250, summary = FALSE)
Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.
## mitchell.resids <- createDHARMa(simulatedResponse = t(preds),
##                             observedResponse = mitchell$Dist_moved/1000,
##                             fittedPredictedResponse = apply(preds, 2, median),
##                             integerResponse = TRUE)
mitchell.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = mitchell_sub$Dist,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
mitchell.resids |> plot()

mitchell.resids |> testDispersion()


    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.49418, p-value = 0.416
alternative hypothesis: two.sided
mitchell.resids <- make_brms_dharma_res(mitchell.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(mitchell.resids)) +
  wrap_elements(~ plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
  wrap_elements(~ plotResiduals(mitchell.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(mitchell.resids))

mitchell.resids |> testTemporalAutocorrelation(time = mitchell_sub$Day)


    Durbin-Watson test

data:  simulationOutput$scaledResiduals ~ 1
DW = 2.0389, p-value = 0.9405
alternative hypothesis: true autocorrelation is not 0
## mitchell.resid1 <- mitchell.resids |>
##   recalculateResiduals(group = mitchell$Day, aggregateBy = mean)
##   ## recalculateResiduals(group = interaction(mitchell$Day,  mitchell$ID),  aggregateBy = mean)
## mitchell.resid1 |> testTemporalAutocorrelation(time=unique(mitchell$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using  as id variables
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
`geom_smooth()` using formula = 'y ~ x'

residuals(mitchell.brm3)[, "Estimate"] |> acf()

Conclusions:

  • the simulated residuals do suggest that there might be a dispersion issue
  • it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.

9 Refit model

## prior(gamma(2, 1), class = 'shape')
priors <- prior(normal(60, 10), class = "Intercept") +
  prior(normal(0, 10), class = "b") +
  ## prior(student_t(3, 0, 10), class = 'sd')  +
  prior(student_t(3, 0, 10), class = "sigma") #+
## prior(lkj_corr_cholesky(1), class = 'cor')
## prior(gamma(0.01, 0.01), class = 'shape')

mitchell.form <- bf(
  Dist ~ scale(Day) + scale(TOD) +
    ar(time = Day, gr = ID, p = 1, cov = TRUE),
  family = gaussian()
)
## family=lognormal())
## family=Gamma(link='log'))
get_prior(mitchell.form, data = mitchell_sub)
                    prior     class     coef group resp dpar nlpar lb ub
                   (flat)        ar                                -1  1
                   (flat)         b                                     
                   (flat)         b scaleDay                            
                   (flat)         b scaleTOD                            
 student_t(3, 61.5, 13.4) Intercept                                     
    student_t(3, 0, 13.4)     sigma                                 0   
       source
      default
      default
 (vectorized)
 (vectorized)
      default
      default
mitchell.brm4 <- brm(mitchell.form,
  data = mitchell_sub,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 1000,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)
Compiling Stan program...
Start sampling
mitchell.brm4 |>
  conditional_effects("Day") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

mitchell.brm4 |>
  conditional_effects("Day") |>
  plot()

mitchell.brm4 |>
  conditional_effects("TOD") |>
  plot()

mitchell.brm4 |>
  conditional_effects("TOD") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

# mitchell.brm4 |> SUYR_prior_and_posterior()

9.0.1 stan plots

The brms package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

  • stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
mitchell.brm4$fit |> stan_trace()

The chains appear well mixed and very similar

  • stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
mitchell.brm4$fit |> stan_ac()

There is no evidence of auto-correlation in the MCMC samples

  • stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
mitchell.brm4$fit |> stan_rhat()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

  • stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

    If the ratios are low, tightening the priors may help.

mitchell.brm4$fit |> stan_ess()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

mitchell.brm4$fit |> stan_dens(separate_chains = TRUE)

9.0.2 pp check

mitchell.brm4 |> pp_check(type = "dens_overlay", ndraws = 100)

9.0.3 DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

  • simulated (predicted) responses associated with each observation.
  • observed values
  • fitted (predicted) responses (averaged) associated with each observation
## mitchell.resids <- make_brms_dharma_res(mitchell.brm4, integerResponse = FALSE)
## wrap_elements(~testUniformity(mitchell.resids)) +
##                wrap_elements(~plotResiduals(mitchell.resids, form = factor(rep(1, nrow(mitchell_sub))))) +
##                wrap_elements(~plotResiduals(mitchell.resids, quantreg = TRUE)) +
##                wrap_elements(~testDispersion(mitchell.resids))

## mitchell.resids |> testTemporalAutocorrelation(time = mitchell_sub$Day)

# mitchell.resids |> testTemporalAutocorrelation(time=unique(mitchell_sub$Day))
autocor_check(mitchell_sub, mitchell.brm3, variable = "Day", n.sim = 250)
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using  as id variables
`geom_smooth()` using formula = 'y ~ x'

autocor_check(mitchell_sub, mitchell.brm4, variable = "Day", n.sim = 250)
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
variog: computing omnidirectional variogram
Using  as id variables
`geom_smooth()` using formula = 'y ~ x'

residuals(mitchell.brm3)[, "Estimate"] |> acf()

residuals(mitchell.brm4)[, "Estimate"] |> acf()

Conclusions:

  • the simulated residuals do suggest that there might be a dispersion issue
  • it might be worth exploring either zero-inflation, a negative binomial model, or including a observation-level random effect.
mitchell.brm3 |>
  augment() |>
  ## group_by(ID) |>
  reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
  ## group_by(ID) |>
  mutate(N = 1:n()) |>
  ## slice(1:12) |>
  ungroup() |>
  group_by(N) |>
  summarise(ACF = mean(ACF)) |>
  ggplot(aes(y = ACF, x = N)) +
  geom_hline(yintercept = 0) +
  geom_segment(aes(yend = 0, xend = N))

mitchell.brm4 |>
  augment() |>
  ## group_by(ID) |>
  reframe(ACF = as.numeric(acf(.resid, plot = FALSE)$acf)) |>
  ## group_by(ID) |>
  mutate(N = 1:n()) |>
  ## slice(1:12) |>
  ungroup() |>
  group_by(N) |>
  summarise(ACF = mean(ACF)) |>
  ggplot(aes(y = ACF, x = N)) +
  geom_hline(yintercept = 0) +
  geom_segment(aes(yend = 0, xend = N))

References

Mitchell, David J, Antoine M Dujon, Christa Beckmann, and Peter A Biro. 2019. Temporal autocorrelation: a neglected factor in the study of behavioral repeatability and plasticity.” Behavioral Ecology 31 (1): 222–31. https://doi.org/10.1093/beheco/arz180.